1. IntroductionIn recent years, high-order harmonic generation (HHG)[1,2] has been a subject of great research interest in strong laser–matter interactions. The HHG has important applications in attosecond science.[4] It can be well described by a three-step model,[5,6] in which harmonics are emitted through the processes of tunneling ionization, propagation, and recombination of the electron with the nuclei. For atoms and symmetric molecules, the emission of harmonics includes only the odd component due to the central symmetry of the atomic and molecular potentials.[7–13] For asymmetric molecules,[14–18] due to the absence of the inversion symmetry, both odd and even harmonics are emitted.[19–24]
Present studies show that odd and even harmonics from asymmetric molecules possess different spectral properties[25–27] and carry different information of the target.[28–33] They therefore need to be studied separately. A simple model[28] has been proposed to describe the generation mechanism of the odd and even harmonics, where the emissions of odd (even) harmonics from asymmetric HeH2+ molecules are considered to be closely associated with the so-called odd (even) transition dipole between the even (odd) component of the asymmetric orbital and the continuum state. This model has been applied to tomographic reconstruction of the asymmetric orbital of the CO molecule with HHG,[33] the analysis of two-center-interference-induced minima in the odd–even HHG spectra,[34] the polarization properties of the odd–even harmonics,[35] the time-resolved study of the odd–even HHG process,[36] and attosecond probing of the vibrational dynamics of asymmetric molecules with odd–even HHG,[37] etc. This model works better for asymmetric molecules with small internuclear distances. For large ones, the resonance between the ground state and the first excited state of the asymmetric system also plays an important role in the HHG process.[38] In this situation, to describe the odd–even HHG, a more complex model which considers this resonance effect is needed.
On the other hand, the HHG spectra from symmetric molecules show a striking minimum,[39] which has been identified as arising from the effect of two-center interference and has attracted broad interest in recent years.[40–44] It has been shown that this minimum in the HHG spectrum corresponds to the minimum in the bound-continuum transition dipole and therefore has important applications in molecular high-order harmonic spectroscopy (HHS). For example, this minimum in the HHG spectrum can be used to read the bond length of the molecule[40,41] and judge the phase of the dipole which is important in the molecular orbital tomography procedure.[45]
For asymmetric molecules, however, in many cases, the striking minimum disappears in the odd or even HHG spectrum and the reason has been attributed to the interplay of different recombination routes originating from the asymmetry of the molecular orbital.[34] Fortunately, with the simple odd–even-HHG model,[28] it is shown that when the minimum cannot be read from the odd–even HHG spectra directly, it can be probed[35] through the polarization measurement[46–53] of the odd–even HHG. Specifically, the position of the minimum in the odd or even dipole corresponds to the harmonic order for the maximal ellipticity of the odd or even harmonics. This polarization measurement of odd–even HHG already shows promise for use in probing the vibrational motion of the asymmetric molecule.[37] However, the prediction of the dipole minimum by the polarization measurement does not always agree well with the theoretical evaluation. In some cases, a remarkable difference is also observed.[35] It is then natural to ask whether this remarkable difference arises from other mechanisms that go beyond the description of the simple model, or if it arises from only the inaccurate evaluation of the dipole in relevant theoretical treatments where some rough approximations are used?
In this paper, we answer the question with improving the evaluation of the dipole. We choose the simplest asymmetric diatomic molecule HeH2+ as the target molecule, with varying molecular parameters such as the effective charge and the internuclear distance. In the frame of linear combination of atomic-orbitals-molecular-orbitals (LCAO-MO) approximation, we use the variational method to obtain the analytical expression for the ground-state wave function of model HeH2+. Using the developed ground-state wave function, we calculate the transition dipole between the ground state and the continuum state. Then we compare the calculated odd–even dipoles with the spectra and ellipticity of odd–even harmonics from the three-dimensional (3D) model HeH2+ obtained through numerical solution of the time-dependent Schrödinger equation (TDSE).
Our simulations show that the use of the variational method to the ground-state wave function remarkably improves the evaluation of the odd–even dipoles, in comparison with the exact ones obtained through diagonalizing the field-free Hamiltonian in one-dimensional (1D) cases. In particular, the minimum in the improved odd–even dipoles agrees well with that predicted by the polarization measurements of odd–even harmonics for different molecular parameters. This agreement also holds as we compare the improved dipoles with the spectra, which are obtained with considering the transition of the continuum electrons back to only the ground state where a striking minimum also emerges. Our results verify the one-to-one matching between the odd–even HHG spectra and the odd–even dipoles, and give important suggestions on the generation mechanism of the odd–even harmonics.
2. Theoretical descriptions2.1. TDSE simulations of odd–even HHGThe Hamiltonian of the model HeH2+ studied here is H(t) = P2/2 + V(r) + r⋅E(t) (in atomic units of ℏ = e = me = 1). We assume that the molecular axis is located in the xoy plane and the laser field is linearly polarized along a direction parallel to the x axis. The potential used here has the form 
. Here 
, 
. Z1 and Z2 are the effective charges. R1 and R2 are the positions of the He and H nuclei to the origin, respectively, with R1 = Z2R/(Z1 + Z2) and R2 = Z1R/(Z1 + Z2). R is the internuclear separation. ξ = 0.5 is the smoothing parameter, and θ denotes the angle between the molecular axis and the laser polarization. For different R, the effective charges Z1 and Z2 are adjusted in such a manner that the ionization potential of model HeH2+ reproduced in our 3D simulations is 1.1 a.u. E(t) = exE(t), where E(t) = f(t)E sin ω0t is the external electric field, and ex is the unit vector along the x axis. f(t) is the envelope function. E and ω0 are the amplitude and the frequency of the external electric field, respectively. In our calculations, we use trapezoidally shaped laser pulses with a total duration of 10 optical cycles and linear ramps of three optical cycles. Numerically, the above Schrödinger equation is solved by the spectral method.[54] We work with a grid of sizes Lx × Ly × Lz = 409.6 a.u. × 51.2 a.u. × 51.2 a.u. for the x, y, and z axes, respectively. The laser wavelength used here is λ = 800 nm and the laser intensity is I = 8 × 1014 W/cm2. In each time step, the TDSE wave function ψ(t) of H(t) is multiplied by a mask function to absorb the continuum wave packet at the boundary. The mask function along the x axis has the form F(x) = cos1/8[π(| x | −x0)/(Lx − 2x0)] for | x |≥ x0 and F(x) = 1 for | x |< x0. Here, x0 is the boundary of the absorbing procedure along the x axis. For the accurate simulations, we have used x0 = Lx/8. The situation is similar for other dimensions of y and z. With the present laser parameters, the maximal classical displacement[5] of the electron along the laser polarization is 
, which can easily be represented in our numerical grids. The absorbing procedure of x0 = Lx/8 is also sufficient to retain the contributions of long–short trajectories and multiple returns[6] to HHG.
Alternatively, we can set the boundary of the absorbing procedure along the x direction (i.e., the direction of the laser polarization) as 
with y0 = Ly/8 and z0 = Lz/8 in the y and z directions unchanged. This treatment removes the contributions of the long trajectory and multiple returns, and the short-trajectory contributions are not influenced. We therefore can analyze the ellipticity of harmonics only arising from the short trajectory. The value of 
corresponds to the maximal displacement of the electron as it travels in the laser field following the short trajectory.
According to the three-step model, it is well known that in each laser cycle, the long and the short electron trajectories contribute significantly to the HHG, which can be modulated through the propagation effect[55] or using phase-stabilized driving pulses.[56] It has been shown that the interference of the long and short electron trajectories[57] greatly influences the ellipticity of the harmonics.[58] As a result, the ellipticity of odd or even harmonics from asymmetric molecules shows strong oscillation and it is difficult to identify the angle dependence of the ellipticity.[35] For these reasons, in this paper, we consider only the short-trajectory TDSE results, for which the elipticity of the odd–even HHG spectra shows a clear dependence on the orientation angle.
Once the TDSE wave function ψ(t) is obtained, the coherent part of the HHG spectrum, parallel or perpendicular to the laser polarization, can be evaluated using
where
ω is the emitted-proton frequency. To study the influence of the excited states on the HHG, the following expression is also used in the evaluation of the HHG spectrum:
which denotes the transition of the continuum electron back to the ground state |0〉 with the amplitude
a0(
t) = 〈0|
ψ(
t)〉. Here, the contributions of the excited states are excluded.
The ellipticity of HHG is determined by the amplitude ratio and the phase difference of the parallel and perpendicular harmonics
where
μ =
S⊥/
S∥ and
δ =
ϕ⊥−
ϕ∥. The intensity and phase of the harmonic components are given by
S∥(⊥) = |
F∥(⊥)(
ω)|
2 and
ϕ∥(⊥)(
ω) = arg[
F∥(⊥)(
ω)]. The range of the ellipticity is 0 ≤
ε ≤ 1. The linear, elliptical, and circular polarizations correspond to
ε = 0, 0 <
ε < 1, and
ε = 1, respectively.
2.2. Analytical description of odd–even HHGAccording to the simple model,[28] the dipole Dodd(ω, θ) that is mainly responsible for the emission of odd harmonics from HeH2+ along the laser polarization e∥ can be written as
with
Godd(
ω,
θ) =
a1 cos(
pkR1 cos
θ) +
a2 cos(
pkR2 cos
θ). Similarly, the dipole
Deven ≡
Deven(
ω,
θ) that is related to the emission of even harmonics along the laser polarization can be written as
with
Geven(
ω,
θ) =
a1 sin(
pkR1 cos
θ) −
a2 sin(
pkR2 cos
θ). Here,
,
, and
.
R1 =
Z2R/(
Z1 +
Z2) and
R2 =
Z1R/(
Z1 +
Z2).
R is the internuclear separation.
Z1 and
Z2 are the effective changes of the He and H nuclei.
Ip is the ground-state ionization potential of HeH
2+.
pk is the effective momentum of the continuum state, |
p〉 ∝ |e
ipk⋅
r〉, with
pk = |
pk| = [2(
Ip +
Ep)]
1/2 that considers the Coulomb acceleration.
[45]Ep is the energy of the continuum state |
p〉, agreeing with the energy conservation relation
Ep =
ω −
Ip. In the expressions of Eqs. (
4) and (
5), the ground-state wave function
ϕ1σ(
r) of HeH
2+ in the mass-center coordinate is roughly approximated by
ϕ1σ(
r) =
Nf[
a1 e
−κra +
a2 e
−κrb]. Here,
Nf is the normalization factor,
ra = |
r −
R1|, and
rb = |
r −
R2|.
In this paper, we improve the description of ϕ1σ(r) with the variational method. According to the LCAO-MO approximation, we assume that the ground-state wave function ϕ1σ(r) of HeH2+ in the mass-center coordinate has the form ϕ1σ(r) = αφa + βφb. Here, α and β are the coefficients. 
and 
are the wave functions of He+ ion and H atom, respectively, with κ1 = Z1 and κ2 = Z2. We approximate the continuum state |p〉 by the plane wave as in Eqs. (4) and (5). Then the dipole Dodd(even)(ω, θ) can be written as
with
D1 ≡ ∫ d
r[e
−κ1re∥⋅
r sin(
pk⋅
r)] and
D2 ≡ ∫d
r[e
−κ2re∥⋅
r sin(
pk⋅
r)]. In the following, we call the dipoles obtained using Eqs. (
4) and (
5) the previous dipoles and call the dipoles obtained using Eqs. (
6) and (
7) the improved dipoles. The calculations of the parameters
α and
β for the ground-state wave function of HeH
2+ with the variational method are presented below.
2.3. Application of variational method to HeH2+The Hamiltonian of the asymmetric molecule HeH2+ is
Here,
V(
r) = −
Z1/
ra −
Z2/
rb is the coulomb potential with
ra and
rb being the positions of the electron to the two nuclei as defined following Eq. (
5).
Z1 and
Z2 are the effective charges, and
R is the internuclear separation. In the frame of LCAO-MO approximation, the bound wave function of HeH
2+ in the mass-center coordinate can be written as
The energy of the asymmetric system can be evaluated using
Substituting Eqs. (
8) and (
9) into Eq. (
10), we have
with
where
Ea = −
Z12/2 and
Eb = −
Z22/2 are the ground-state energies of He
1 ion and H atom, respectively. According to the normalization condition for
φa and
φb, we have
Saa =
Sbb = 1 and
Sab =
Sba. Taking the derivative of
α or
β in Eq. (
11), we obtain the following equations:
The secular equation for Eq. (
12) is
To solve the secular equation, we have
AE2 +
BE +
C = 0, with
,
B = −(
Haa +
Hbb) + (
Hab +
Hba)
Sab, and
. The solutions of the above equation are the energies of the 1
σ and 2
σ states
According to the normalization condition of Eq. (
9), we have
Then we arrive at
Substituting Eqs. (
12) and (
13) (Eq. (
14)) into Eq. (
15), we have
where
Then we obtain the 1
σ and 2
σ wave functions of HeH
2+
3. Results and discussionIn the following, we apply the variational method to model HeH2+ with different molecular parameters such as the effective charges Z1 and Z2 and the internuclear distance R. We make comparisons between the previous dipoles and the improved dipoles. Then we use these dipoles as a benchmark to analyze the HHG mechanism of odd–even harmonics from asymmetric molecules.
3.1. 1σ-state and 2σ-state energy curves and improved ground-state wave functionIn Fig. 1, we compare the R-dependent ground-state energy of HeH2+ obtained by the variational method (black-triangle) with that obtained by the exact numerical simulations (red-circle). The results of the two different methods show good agreement with each other.
In Table 1, we show the coefficients of α and β obtained with the variational method for HeH2+ at different R. As shown in Table 1, the parameter α1 is much larger than β1, implying that for the 1σ state of HeH2+, the probability density distribution around the helium nucleus is much higher than that around the hydrogen nucleus. In contrast, for the 2σ state, the electron is mainly located around the H nucleus (the parameter α2 is smaller than β2).
Table 1.
Table 1.
 Table 1.
Coefficients of the wave functions of Eq. (18) for HeH2+ obtained using the variational method with Z1 = 2 and Z2 = 1 at different internuclear separations R. α1 and β1 obtained using Eq. (16) are the coefficients for the 1σ state. α2 and β2 obtained using Eq. (17) are the coefficients for the 2σ state. .
| R/a.u. |
α1 |
β1 |
α2 |
β2 |
| 1.0 |
0.9620 |
0.0586 |
−0.8580 |
1.2876 |
| 1.2 |
0.9552 |
0.0756 |
−0.7486 |
1.2112 |
| 1.4 |
0.9586 |
0.0777 |
−0.6481 |
1.1545 |
| 1.6 |
0.9659 |
0.0723 |
−0.5572 |
1.1128 |
| 1.8 |
0.9738 |
0.0631 |
−0.4762 |
1.0821 |
| 2.0 |
0.9810 |
0.0531 |
−0.4043 |
1.0597 |
| 2.2 |
0.9867 |
0.0434 |
−0.3417 |
1.0433 |
| 2.4 |
0.9909 |
0.0348 |
−0.2878 |
1.0312 |
| 2.6 |
0.9929 |
0.0274 |
−0.2416 |
1.0224 |
| 2.8 |
0.9960 |
0.0215 |
−0.2023 |
1.0161 |
| 3.0 |
0.9973 |
0.0167 |
−0.1691 |
1.0114 |
| 3.2 |
0.9983 |
0.0129 |
−0.1410 |
1.0081 |
| 3.4 |
0.9989 |
0.0099 |
−0.1175 |
1.0057 |
| 3.6 |
0.9993 |
0.0077 |
−0.0978 |
1.0040 |
| 3.8 |
0.9995 |
0.0059 |
−0.0812 |
1.0028 |
| 4.0 |
0.9996 |
0.0045 |
−0.0674 |
1.0020 |
| 4.2 |
0.9998 |
0.0035 |
−0.0559 |
1.0013 |
| 4.4 |
0.9999 |
0.0026 |
−0.0463 |
1.0009 |
| 4.6 |
1.0000 |
0.0020 |
−0.0383 |
1.0007 |
| 4.8 |
1.0000 |
0.0016 |
−0.0317 |
1.0005 |
| 5.0 |
1.0000 |
0.0012 |
−0.0262 |
1.0003 |
| Table 1.
Coefficients of the wave functions of Eq. (18) for HeH2+ obtained using the variational method with Z1 = 2 and Z2 = 1 at different internuclear separations R. α1 and β1 obtained using Eq. (16) are the coefficients for the 1σ state. α2 and β2 obtained using Eq. (17) are the coefficients for the 2σ state. . |
These coefficients also show that the 1σ wave function ϕ1σ(r) = αφa + βφb obtained with the variational method differs remarkably from the previous one ϕ1σ(r) = Nf[a1 e−κra + a2 e−κrb. One can expect that this difference will play an important role in the calculation of the odd–even dipoles, as shown below.
The accurate calculation of the bound-continuum transition dipole 〈0|r|p〉 in 3D cases is difficult, as it needs the knowledge of both the bound state |0〉 and the continuum state |p〉. Here, we evaluate the exact dipole in 1D cases where all of the bound and continuum states of the system can be obtained by diagonalizing the field-free Hamiltonian 
at the fixed-nuclei approximation. The 1D asymmetric Coulomb potential used here has the form 
, with ξ = 0.5 and Z2 = Z1/2. R1 = Z2R/(Z1 + Z2) and R2 = Z1R/(Z1 + Z2). For different internuclear distances R, we adjust the parameter Z1 in such a manner that the ionization potential of the 1D asymmetric system reproduced in our numerical simulations equals 2.25 a.u. (the ionization potential of HeH
). When the calculated continuum state has an odd-like (even-like) parity, through the expression 〈0|r|p〉, we obtain the corresponding odd (even) exact dipole,[28] as shown in Figs. 2(a) and 2(b). The parameters Z1 and Z2 used in the numerical simulations are also used in evaluating the analytical dipoles with Eqs. (4)–(7) in Fig. 2.
One can observe from the intersections of the odd (solid-black) versus even (dashed-red) dipoles in Fig. 2 that the improved dipoles agree better with the exact ones. We use the vertical-dashed arrows to indicate several of these intersections. For the positions of the minima in the odd or even dipoles, the improved ones are also closer to the exact ones, especially for the minima in the odd dipoles, as indicated by the vertical-solid arrows in the right column of Fig. 2. These results verify the applicability of the ground-state wave function obtained with the variational method in calculating the odd–even dipoles for HeH2+.
It should be stressed that although the use of the variational method improves the agreement between the exact and the analytical odd–even dipoles. There are still some differences between them, both in the positions of the minima in the dipoles and the intersections of the odd versus even dipoles. These differences can arise from the plane-wave approximation with the effective momentum pk[59] for the continuum electron in calculating the dipoles, where the Coulomb effect is not well described. One can expect that this Coulomb effect is stronger for molecules with larger ionization potentials.
3.2. Comparison between dipoles and spectraIt has been shown[35] that for model HeH2+ with small internuclear distances R, the previous dipole of Eq. (4) associated with odd harmonics shows a striking minimum. However, the minimum is absent in the odd TDSE spectrum. By contrast, the ellipticity curve of odd harmonics shows a maximum at one harmonic order at which the minimum appears in the odd dipole. The ellipticity of harmonics therefore can be used as a tool to probe the position of the minimum in the dipole. This dipole minimum is important, as it encodes the information of the molecular structure and has potential applications in asymmetric molecular orbital imaging.[33] However, in some cases, the position of the minimum in the dipole does not agree well with the harmonic order for maximal ellipticity. To use the ellipticity measurement of harmonics as a tool to judge the dipole minimum, it is necessary to clarify this disagreement.
In Fig. 3, we compare the spectra, ellipticity, and dipoles of odd harmonics for a model HeH2+ molecule with Ip = 1.1 a.u., Z1/Z2 = 1.5, and R = 1.5 a.u. at different orientation angles θ. The spectra are obtained through 3D TDSE simulations.
First, one can observe that the accurate TDSE spectra of Eq. (1) (solid-black) in the first row of Fig. 3 do not show a striking minimum. By contrast, the spectra approximated with Eq. (2) (dashed-red), where the transition of the continuum electron back to only the ground state is considered, show a striking minimum. The position of the spectral minimum shifts towards higher harmonic orders as the orientation angle increases. The comparisons imply that the transition of the continuum electron back to the excited states plays an important role in the HHG of asymmetric molecules, resulting in the disappearance of the minimum in the accurate TDSE spectrum of Eq. (1).
Secondly, corresponding to the spectral minima in the first row of Fig. 3, the ellipticity curves of harmonics in the second row of Fig. 3 show a striking maximum. In particular, the maximal ellipticity appears at the harmonic order at which the improved dipoles of Eq. (6) (dashed-red) in the third row of Fig. 3 show a striking minimum, as indicated by the vertical arrows. In contrast, the positions of the minima in the previous dipoles of Eq. (4) (solid-black) in the third row of Fig. 3 differ remarkably from the improved ones and this difference is more remarkable for larger orientation angles θ. These comparisons reveal that there is a one-to-one matching between the dipole minimum and the ellipticity maximum. The previous dipoles of Eq. (4) with a rough approximation for the ground-state wave function underestimate the positions of the minima for about 10 to 30 harmonic orders here.
To check our results, we also perform simulations at other molecular parameters, as shown in Figs. 4 and 5, where we fix the value of Z1/Z2 and change the value of the internuclear distance R. On the whole, in all cases, the position of the minimum in the improved dipole agrees with the harmonic order at which the maximal ellipticity emerges. By contrast, the position of the minimum in the pervious dipoles differs remarkably from the harmonic order of the maximal ellipticity in some cases. It should be mentioned that for the case of the large angle of θ = 50° in Fig. 4(i), the improved dipole has no minimum. The corresponding HHG spectra in Fig. 4(c) also do not show a minimum. In this case, the correspondence between the improved dipole and the ellipticity is not very obvious. However, as we increase the laser wavelength (such as λ = 900 nm) with extending HHG cutoff, both minima appear in the HHG spectra and the improved dipole, and the correspondence between the improved dipole and the ellipticity becomes striking. In addition, we also fix the value of R and change the value of Z1/Z2, the close relation between the dipoles, spectra, and ellipticity is also observed in our extended simulations. For simplicity, in the above discussion, we have chosen the HeH2+ with an active electron as the target molecule. This close relation is also expected to appear for other asymmetric molecules such as CO with more electrons, for which the odd–even dipoles also give a good prediction of the relative yields of odd–even harmonics.[33] When applying the variational method to CO to obtain the analytical expression of the 1σ ground-state wave function, the multiple-electron effect needs to be considered.
